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Abstract 

We introduce three models of fragmentation in which the largest fragment in the system can 
be broken at each time step with a fixed probability, p. We solve these models exactly in the 
long time limit to reveal stable time invariant (scaling) solutions which depend on p and the 
precise details of the fragmentation process. Various features of these models are compared 
with those of conventional fragmentation models 
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1 Introduction 



The generation of random fragments by successive and sequential breaking is an irreversible 
kinetic process that occurs in a wide range of physical systems in both science and technology 
[1-6]. This variety of apphcations has motivated theorists to obtain analytical solutions to 

simple fragmentation models as a function of both the fragmentation rules and the initial 
conditions [7-10]. 

In this letter we introduce three unusual models of fragmentation which, at each time step, 
allow the fragmentation of the largest fragment in the system with some externally tuneable 
probability. In Model A the largest fragment is broken into two pieces at every time step. 
The size of the resultant fragments is determined by an arbitrary probability distribution. 
In Model B, the largest fragment is split with probability p and with probability (1 — p) 
another fragment is chosen. When the later choice is made, every fragment is equally like to 
be chosen. In Model C, this choice is controlled by a homogeneity index /3. These models are 
solved exactly in the long time limit to reveal stable time invariant probability distributions. 
The results are compared with conventional models of fragmentation. 

This work was originally motivated by three different observations. Firstly, the extremal 
properties of fragmenting systems have been the subject of interest recently [11,12]. This 
work considered systems in which a fragment is put to one side at each time step and takes 
no further part in the fragmentation processes. This process leads to singularities in the 
probability distribution for the largest fragment [11]. 

Secondly, the kinetics of systems of growing and coalescing (equal-sized) droplets has 
been considered in which the two droplets which were closest together coalesced at the next 
time step [13]. This lead to some interesting kinetic behaviour, in particular the droplet 
patterns that emerged were self-similar in time with inter-droplet distances obeying simple 
scaling behaviour. We wondered, if the analogous system in fragmentation models would 
yield similarly interesting results. 

Finally, we observed that the largest particle was always the most likely to fragment in 
any non-shattering fragmentation process. If we allowed the largest fragment to be chosen 
for fragmentation with some fixed probability this would allow us to utilise the methods 
in [13] to obtain solutions to the kinetic process. This fixed probability may seem a little 
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unrealistic, but it is not difficult to envisage a kinetic system, such as froth under some 
external pressure, in which the largest bubble, with its greater volume and surface area, is 
much more likely to fragment than any of the other bubbles in the system. 

In the remainder of this introduction we give a brief summary of previous work on models 
of binary fragmentation. Then in the next three sections we introduce the models and present 
our analytical solutions. In section 5, we discuss the scaling behaviour of the models and in 
the final section summarise our findings and conclusions. 

In most of these models of binary fragmentation, the concentration of fragments of size 
y at time t, n{y,t), evolves according to 



n{y,t + 6t) = n{y,t) — 6tn{y,t) / R{z,y — z)dz + 26t) / R{y, z — y)n{z,t)dz (1) 



where, R{y, z) is the intrinsic rate at which a particle of size {y + z) breaks into particles of 
size y and z. The first term on the right hand side is the contribution from those particles not 
chosen for fragmentation in time {t,t + 6t). The second term represents the decrease in the 
number of particles of size y by fragmentation into particles of size z and {y — z){y > z). The 
third term represents the increase in the number of particles of size y due to the fragmentation 
of particles of size z{> y), such that one of the products is of size y. 

The upper limit of the second integral on the right hand side is usually set to some fixed 
value greater than the size of the largest particle in the system. This value is usually 1 or 
oo depending on the details of the model. To connect with the analysis of our models in the 
following sections, we have set the upper limit to L{t), the size of the largest particle at time 
t. As n{y, t) = for all y > L{t), this choice has no effect on our results. 

A general scaling theory has been constructed for systems of fragmentation particles [10]; 
the particle size distribution in the limit y ^ 0,t ^ oo can be written as 



where s{t) is the time dependent typical cluster mass and the exponent 2 is required by mass 
conservation. Exact solutions have been obtained to the kinetic equation (1) with a number 
of different fragmentation rules R{y,z). In particular, solutions have been obtained [3,8,9] 
when R{y, z) — {y + z)^~^. This describes a situation in which the rate of fragmentation is 
determined solely by the size of the particle and every point along the line of the particle is 





(2) 
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equally likely to fragment. When a; > this system exhibits simple scaling dynamic scaling 
and s{t) ~ with z — For a; < shattering occurs and a cascading break-up of 

particles causes a finite amount of mass to be transferred to an infinite number of particles 
with zero mass [8,10]. We now turn to detailed consideration of our models. 

2 Model A 

Firstly, we consider a model in which the largest particle is fragmented into two pieces at 
each time step. If we assume that in time (t, t + 5t) the size of the largest fragment changes 
from L{t) to {L{t) — 5L) then we can choose R{y, z) as 

5tR{y, z) = ^^fi^^)Siy + z - m) (3) 

This term is easily understood; the (^-function ensures that only the largest particle is frag- 
mented, f{x) is defined as the rate with which particles of size xL and {1 — x)L are created 
from the largest particle L{t) and is the probability of placing the cut in the largest 
particle in a particular infinitesimal length 6L. This rate is time dependent, through its 
dependence on L{t), which is in line with the definition of the model. Inserting the rate into 
the kinetic equation (1) yields 

niy, t + 5t)^ n{y, t) + ^n{L{t),t)f{^) (4) 

where the second term on the right-hand side describes the gain in the number of particles 
of size y from the fragmentation of the largest fragment. 

We can define the normahsed density of fragments of length y at time t as 

g(^y t) = ^(^'^) 

Jo^^^ n{z,t)dz 

so that, as L{t) — > L{t) — 6L in time t ^ t + 6t, 



niy,tSt) 
jQ^''~^^n{z,t + dt)dz 



9(v.t+H)^ .... ■ w 



Consequently, following the methods of [13], we can rewrite equation (4) as 

g{y, t+5t)^ g{y, t)[l - 5L^(L, t)] + ^^(L, (7) 
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The natural length scale in this model is L{t), the length of the largest fragment. So, to 
solve equation (7), we take account of this by introducing F{x,t), defined by, 

F{x,t)^L{t)g{xL{t),t) (8) 

Notice that x, defined hy x — y/L{t), is restricted to the range [0, 1]. Remembering that 
L{t) is a function of time , we can obtain a partial differential equation for F{x, t); 

= 2F(1, t)f(x) - F(l, t)F(x, t) - F(x, t) - x^^^ (9) 

Here, we have set Sx/St, which determines the relationship between the real time and the 
length of the largest fragment L{t), to 1. This choice is arbitrary and has no effect on the 
kinetics of the process. 

The remainder of this section is concerned with solutions to this equation. First of all 
we will assume that F{x,t) evolves to a time independent quantity as t ^ oo. In this limit 
we can solve (9) to yield 

Fix) = ^ - 2^ £ my'^'^dy (10) 

where F{x) — limj^oo F{x, t). 

There are a number of constraints on f{x). Obviously to preserve symmetry we must 
have f{x) — f{l — x) and the distribution is normalised so that 



ff{y)dy^l. (11) 
Jo 



These two constraints lead to 



[\f{y)dy^ 1/2. (12) 
Jo 

A choice of f{x) which satisfies both these constraints is 

Jo ?/"(! - yrdy 

where wc assume a > —3/2 so that the variance of f{x) is well defined. Inserting (13) into 
(10) and integrating F{x) between and 1 reveals 

t F{x)dx = 1 + [1 - ^P±±Imll!l±}l^^ lin, .-(1) (14) 
Jo /3{a + l,a + l) ^x-^o ^ ' 



4 



For F{x) to be correctly normalised the second term on the right hand side must be zero. 
The condition F(l) > 1 means that this term diverges unless the contents of the square 
bracket are zero. It is a simple matter to show that this can be achieved for all a by setting 
F(l) = 1. As a result, we rewrite equation (10) as 

F{x)^^[l-2 [\f{y)dy]. (15) 

Substituting (13) into equation (15) and evaluating the integral reveals that for a = we 
have the trivial solution F{x) = 1 for all x. For a = 1 we have 

F(x) ^x{A-3x). (16) 

for a — 2, 

F(x) = x^(15 - 2Ax + lOx^) (17) 

and for a = —1/2 

F{x) = 4 - -^[Cos-\V^) + ^x{l-x)]. (18) 

We checked these solutions numerically by calculating F{x, to) at some (long) time to for a 
number of realisations of the randomness in the system and then averaging F[x,to) over a 
large number of these realisations. The average values of F{x, to) we obtained agreed with 
the analytical solutions in equations (15)-(18). 

We can also consider the behaviour as x — > 0. If we assume that X ^ 0,f{x) ^ Ax^ 
then we find that 

9 A 

n^) - (19) 

in the same limit. If we take f{x) to have the form (13) then we have /3 — a and A — 
r{2a + 2)/ (r(a))^. Alos, notice that when a < 0, F{x) diverges as a; ^ and when a > 0, 
F{x) goes to zero smoothly. This behaviour in the probability distribution at x = is similar 
to that observed in other models [3,10]. 

The above analysis is restricted to the stationary state. However, as these models are rel- 
atively simple, one can attempt to obtain explicit time dependent solutions. For instance, us- 
ing the method of characteristics, one can solve equation (9) when f{x) = 1 self-consistently 
to reveal 

Fix, t) = 2 F(l, i^)e"-*H + F(xe-\ 0)e-*+Jo (20) 

JQ 



One can then in principle, insert some initial conditions F{x,0), to solve for f{l,t) then 
substitute F{x, 0) and F(l, t) into (20) to give the exphcit solution. The only case we have 
been able to do this when F{x,0) = 1 for < x < 1. This yields the trivial stationary 
solution F{x,t) = 1. We were not able to solve (9) by the method of characteristics for any 
more complicated f{x). 

3 Model B 

In this section we consider a model in which, at each time step, either the largest fragment 
is selected with probability p or another fragment is chosen with probability (1 — p). 

If the largest fragment is chosen, it is split into pieces xL and {l — x)L with a probability 
independent of x. That is to say the daughter distribution function is uniform, or in the 
formalism of Model A, f {x) — 1 or a — 0. If it is decided to split a smaller fragment, every 
fragment is equally likely to be chosen for fragmentation. This corresponds, to the onset of 
shattering in conventional fragmentation models when a; = 0. Consequently, when p — 1, 
this model is equivalent to Model A, when p = 0, it is equivalent to a model of fragmentation 
originally solved by Ziff and McGrady [9]. 

The fragmentation rate for this model is given by 

StR{y, z) = {l-p)^+ p^d{y + z - L{t)) (21) 
y + z L{t) 

and the kinetic equation for the process is 

6L /-^W dz 

n{y, t + 5t)^ niy, t) + 2p-—n{L{t),t) - (1 - p)5tn{y, t) + 2(1 - p)5t / n{z, t) — . (22) 

L{t) Jy Z 

Using an identical method to that used for Model A we can obtain 

= 2F(1, t) - [3 - 2p + F(l, t)]F{x, t) - x^^^ + 2(1 - p) £ F{y, t)^ (23) 

where we have set p5x/5t = 1 and assumed that p > 0. In the limit i ^ oo we can solve 
this equation self-consistently for p > to reveal the time invariant solution 

F{x) = px^-^ (24) 

Here, in contrast to Model A, F(l) is picked out by the solution of the equation and there is 
no need to consider the normahsation of F{x) to find -F(l). As before, this result has been 
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confirmed by numerical simulation. Models A and B coincide with the same special, trivial 
solution when a — and p — 1. 



4 Model C 

In this section we consider a generalised version of Model B in which the rate with which 
the non-largest fragments are chosen for fragmentation is controlled by a homogeneity index 
(3. More precisely, a particle of size y + z is split into fragments of size y and z with a rate 
(y + z)^~^ [9]. This equivalent to cu — (3 in conventional models. Consequently, R{y,z) is 
given by 

5tR{y, z) = {l- p)5t{y + zf-^ + P^^V + z - L{t)) (25) 
and the rate equation for this process is 

n{y, t + St) = n{y, t) + 2pj^n{L{t),t)) - (1 - p)5ty^n{y, t) + 2(1 - p)5t n{z, t)z^-^dz 

(26) 

In this model, if one sets /3 = 0, we recover Model B and when p = 1, we get Model A. 
We now follow the same steps as before to reveal an integro-differential equation for 

F{x,t)- 

= 27F(1, t) - [7(1 + F(l, t)) + (1 - p){xP + /' F(y, i)/dy)]F(x, t) 

Jo 

- 7^^^^ + 2(1 - p) [ F{y, t)y^+'dy (27) 



dt 



where the time dependent function 7(t) is defined by 

P^-^ = L{tr^{t). (28) 

We choose the relationship between the real time t and the largest fragment L{t) so that in 
the limit t oo, j{t) goes to some constant, to be determined, 7. Consequently, in the long 
time limit when the partial time derivative of F{x,t) goes to zero, equation (27) becomes 

^x^^ + [^{l+F{l)) + {l-p){x^ + £ F{y)y^dy)]F{x) = 27F(l) + 2(l-p) jy{y)y^-^'dy 

(29) 

This can be solved self-consistently using the functional form for F[x)] 

F(x) ^ ae-'"'" . (30) 



Substituting (30) into (29) reveals two equations for a and h, 

J^e-Hi^-^dt 

and 



(31) 



where the integral in the denominator of a is an incomplete gamma function. As in Model 
B, we find that F{1) — p is picked out as a condition for a function of the form (30) to be 
a solution. Equations (31) and (32), together with F(l) = p, completely determine a, b and 
7. So, as an example, we consider P — 1, when we obtain an exponential F{x) 

F{x) = pe''(^-^) (33) 

where 

b^p(e''-l). (34) 

Wc have confirmed this solution, and those for some other values of /?, by numerical simu- 
lation. 

5 Scaling 

In this section we consider the scahng properties of all three models. The function F{x,t) 
is related to the number density of the fragments, n{y, t), via 

where, N{t) is the total number of fragments. Notice that the total length of the fragments, 
a constant, is given by 

/ yn(y,t)dy ^ N(t)L(t) xF(x,t)dx. (36) 
Jo Jo 

In the limit t ^ 00, when the integral on the right hand side becomes time independent, we 
know that N{t) ~ t, so we can deduce that L{t) ~ 1/i. Inserting this result into equation 
(35) gives the scaling form, 

n{y,t)r.t'F{^)^f^yt) (37) 



indicating that the typical cluster mass s{t) ~ 1/t and the kinetic exponent 2; = 1 for all 
these models. So, in the scaling hmit i — > 00, the scaling function $(^) ~ F{^), where 

C ~ yt- 

To illustrate this scaling result, and to show the agreement with our analytical result, we 
have done two numerical simulations. In figure 1, we have shown a plot of n{y, t)/t^ against 
yt from a simulation of Model B with p = 1/3 for three different times. This figure can be 
thought of as a plot of $(^) against ^ once each axis has been rescaled by constant prefactors. 
We see that the data collapses as predicted and in line with our analytical results (equation 
(19)), $(0 ^ for 'C < 1- The analytical result is shown as solid line on the figure. The 
region ^ > 1, where $ is zero and the point ^ = 1, where $ goes to zero discontinuously, is 
not shown. 

In figure 2 there is similar plot for Model A with a = 1. After rescaling both axes 
by suitable prefactors this result exactly matches the analytic prediction in equation (13); 
^(0 ~ C(4 — 3^) for ^ < 1 and $(^) = for ^ > 1. This prediction is shown with the solid 
line. In this figure we have included the region around y ^ L{t) ^ 1/t (about yt — 0.034 
in the units on the figure, or ^ = 1) where the data collapse is not quite so good. This is 
due to the difficulty of capturing the discontinuity at y ^ L{t){^ = 1) numerically. However, 
one can see that the slope of the curve around y ^ L{t) gets progressively steeper as t gets 
larger. This is as expected; the discontinuity would be captured exactly in the limit i — > 00. 

In the above discussion we saw that the largest fragment L{t) ~ t"^ with 5 — 1 in the 
long time limit. To compare this particular result with conventional models of fragmentation 
in which a particle of size y + z is split into fragments of size y and z with a rate equal to 
1/ {y -\- z){uj — 0) or independent of y and z{lo = \). In these two models, we found that the 
size of the largest fragments also decays algebraically, but more slowly, as 5 with 5 = 0.20 
and 0.85 for a; = and 1 respectively. 

6 Summary and Conclusions 

To briefly summarise our results, we have introduced three models of fragmentation in which 
there is a fixed probability per time step of breaking the largest fragment in the system. In 
Model A the largest fragment was fragmented at each time step with an arbitrary fragmen- 
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tation kernel. In Model B and C the largest fragment was chosen for fragmentation with 
probability p and (1 — p) another fragment was chosen. The behaviour of these models is 
completely described by the function F{x,t), which evolves to a stable distribution F{x) in 
the limit t ^ oo. We were able to obtain an analytic solution for F[x) for a number of 
different fragmentation kernels. We found that the models exhibited simple scaling , with 
kinetic exponent z — 1 and the scaling function was equal to F{x). 

To study the kinetics of these models we looked for solutions of the kinetic equation 
obeyed by F{x,t). In the simplest case of Model A (equation (9)), this is a non-linear in 
t, non-local in x, partial equation of the two variables t and x. This means that solutions 
must be obtained self-consistently and this limits the scope for finding exact time dependent 
solutions to all but the most trivial problems. Our inability to study the kinetics analytically 
leads one to wonder about the range of validity of the stationary solutions we have obtained. 
Do all initial conditions lead to these stationary solutions, or if not, how is the phase space 
of initial conditions divided up? We have performed a number of numerical simulations 
and have always found that these models evolve to the predicted stationary states. We 
suspect that all well behaved, physically important, smooth initial conditions evolve to these 
stationary states. However, if we are able to find a set of initial conditions that do not evolve 
to these stationary states, we will return to this subject in a future publication. 
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Figure Captions 

1. n{y, t)/t^ against yt for Model B with p— 1/3 for three different times. The sohd 
hne is a plot of the analytical result. 

2. As figure 1 for Model A with a = 1. 
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